Theoretical and numerical investigation of the shock formation of dust ion acoustic 

waves* 
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We present a theoretical and numerical study of the self-steepening and shock formation of large- 
amplitude dust ion-acoustic waves (DIAWs) in dusty plasmas. We compare the non-dispersive 
two fluid model, which predicts the formation of large amplitude compressive and rarefactive dust 
ion-acoustic (DIA) shocks, with Vlasov/fluid simulations where ions are treated kinetically while a 
Boltzmann distribution is assumed for the electrons. 



Shukla and Silin [l| predicted the existence of small 
amplitude dust ion-acoustic waves (DIAWs) in an un- 
magnetized dusty plasma. In the DIAWs, the restor- 
ing force comes from the pressure of inertialess elec- 
trons, while the ion mass provides the inertia to sup- 
port the waves. On the timescale of the DIAWs, charged 
dust grains remain immobile, and they affect the over- 
all quasi-neutrality of the plasma. When the dust grains 
are charged negatively, one has the depletion of the elec- 
trons in the background plasma. Subsequently, the phase 
speed [u/k = (n ia /n e0 + 3Ti/T e ) 1/2 C s , where n ia (n eQ ) is 
the unperturbed ion (electron) number density, T, (T e ) 
is the ion (electron) temperature, C s = (Te/rrii) 1 / 2 is the 
ion acoustic speed, and m, is the ion mass] of the DI- 
AWs becomes larger than the usual ion-acoustic speed in 
an electron-ion plasma without negatively charged dust, 
since n i0 > n eQ . When n iQ ^> n e0 or Tj 3> T e , small- 
amplitude DIAWs do not suffer Landau damping in a 
plasma, since the increased phase speed is much larger 
than the ion thermal speed (Ti/nii) 1 ^ 2 . Small-amplitude 
DIAWs have been observed in laboratory experiments [2| , 
and the observed phase speed is in an excellent agreement 
with the theoretical prediction of Ref. . 

Recently, laboratory experiments 0, 0, 0, have been 
conducted to study the formation of dust ion-acoustic 
(DIA) shocks in dusty plasmas. Dust ion acoustic com- 
pressional pulses have been observed to steepen as they 
travel through a plasma containing negatively charged 
dust grains. Theoretical models 0,0 have been proposed 
to explain the formation of small amplitude DIA shocks 
in terms of the Korteweg-de Vries-Burgers equation, in 
which the dissipative terms comes from the dust charge 
perturbations j^] ■ Popel et al. ^(J have included sources 
and sinks in the ion continuity equation, linear ion pres- 
sure gradients in the nonlinear ion momentum equation 
with a model collision term, as well as the dust grain 
charging equation to study the formation DIA shock- 
wave structures. 
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In this Brief Communication, we present analytical and 
numerical studies of large amplitude DIA shock waves 
in an unmagnetized dusty plasma [llj. We use fully 
nonlinear continuity and momentum equations for the 
warm ion fluid, as well as Boltzmann distributed elec- 
trons and the quasi-neutrality condition to examine the 
spatio-temporal evolution of large amplitude dust ion- 
acoustic pulses. We find simple-wave solutions of our 
fully nonlinear two fluid model, and compare them with 
those deduced from the time-dependent Vlasov simula- 
tions which uses initial conditions corresponding to the 
ones obtained from our theoretical model. 

We consider an unmagnetized dusty plasma whose 
constituents are singly charged positive ions, electrons 
and charged dust grains. Thus, at equilibrium, we have 
Uio = n e Q — eZdndo, where e equals —1 (+1) for negatively 
(positively) charged dust grains, Z d is the number of el- 
ementary charges residing on the dust grain, and n^o is 
the equilibrium dust number density. On the timescale of 
our interest, the dust grains are assumed to be immobile. 
The dynamics of low phase speed (in comparison with 
the electron thermal speed) nonlinear, dust ion-acoustic 
waves is governed by a Boltzmann distribution for the 
electrons 
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where n e (m) is the total electron (ion) number density, 
e is the magnitude of the electron charge, <f> is the wave 
potential, and Vi is the ion fluid velocity. The system is 
closed by means of Poisson's equation 
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In the following, we consider non-dispersive DIAWs, 
and use the quasi-neutrality condition n e = + eZ^n^o 
instead of Eq. (4), together with the normalized variables 
N = ni/riio, u = Vi/C s , ip = ecj>/T e , z = r D x x and 
t = ujpit, where uj p i = (A-Kmoe 2 jrai) 1 ! 2 is the ion plasma 
frequency and rp — C s /uj p i is the electron Debye radius. 
Thus, the system of equations (l)-(3) can be rewritten as 
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where a — n^jriio and r\ — Tio/T e . In obtaining Eq. (5), 
we have used ip = ln[(iV + a — l)/a] which follows from 
Eq. (1) and the quasineutrality condition. 

In order to study the nonlinear evolution of large am- 
plitude DIAWs, we seek simple wave solutions [IJ of 
Eqs. (5) and (6). For this purpose, we rewrite them in 
the matrix form as 
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Here, the nonlinear wave speeds are given by the eigen- 
values 



A± = u±N 1/2 
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of the square matrix multiplying the second term in Eq. 
(7). The square matrix in Eq. (7), which we denote A, 
can be diagonalized by a diagonalizing matrix C whose 
columns are the eigenvectors of A, so that 




FIG. 1: The wave speed A+ as a function of N for different 
values on r\ and a. 
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where the new variables are 
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Equations (12) and (13) describe the DIAWs propagat- 
ing in the positive and negative z directions, respectively. 
Setting ?/>_ to zero, we have 
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Multiplying Eq. (7) by C" 1 from the left gives the 
diagonalized system of equations 



dr 



dz 



o. 



(12) 



u(N) 



N 



[l/N'(N' + a- l)+3r)] 1/2 dN', (15) 



which inserted into Eq. (12) gives 

A +< w » - f (swr^)"')"' * N ' +Nl '* («rbr 

(16) 

Since can be written as a function of N, Eq. (12) 
holds also for N, i.e. 



(17) 



which, as long as N is continuous, has the general so- 
lution N = /o(£), where — x — X + (N)t and fo is the 
initial condition for N. We have plotted A+ as a func- 
tion of N in Fig. 1. Here, we see that A + grows with 
increasing N in the two cases with a = 0.25. In the case 
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with A = 0.05, however, the phase speed first decreases 
for N 1 before increasing with increasing N. In the 
small-amplitude limit, viz. N — l+Ni, where |JVi| << 1, 
we have the first-order Taylor expansion A + = c + jNi , 
where c = (l/a + 377) 1 / 2 is the linear acoustic speed and 
7 = (3a+ 12rja 2 — l)/2a(a + 3ria 2 ) 1 / 2 is the coefficient in 
front of the nonlinear term. We note that 7 is negative 
for sufficiently small a (in agreement with the a = 0.05 
case displayed in Fig. 1), and in the cold ion limit (r/ = 0) 
we recover the result that a < 1 /3 leads to a negative co- 
efficient 0,0] in front of the nonlinear term. The linear 
acoustic speed increases when a decreases. Thus, in the 
presence of negatively charged dust, the phase speed of 
the waves may becomes much larger than the ion acous- 
tic speed, so that the Landau damping of the waves de- 
creases pj. 

In order to compare the fluid and kinetic theories, we 
have solved the coupled Eqs. (5) and (6) numerically 
and compared the results with numerical solutions of the 
Vlasov equation. As an initial condition for our fluid 
simulations, we take a large-amplitude localized density 
pulse, N = 1.5 - 0.5sech[3sin(27rz/20000) + 1.5], while 
the initial condition for the velocity is obtained from the 
simple wave solution as u(N) = [1/N'(N' + a — 1) + 
3rj] x l 2 dN'. The results are compared with numerical so- 
lutions of the ion Vlasov equation 



dr dz dz dv 



(18) 



where v has been normalized by C s and the ion distribu- 
tion function / by n^o/Cg. Here, we have also used the 
quasineutrality condition and thus ip — \n{(N+a — I) /a], 
where N = f fdv. For the initial condition, we are 
using the shifted Maxwellian ion distribution function 
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where we are using the same initial condition for the den- 
sity N and the velocity u as in the fluid simulations. For 
the scaled (by X^o) ion temperature T, we obtain an ini- 
tial condition by combining the ideal gas law P = NT, 
where P = Pi/Pio is the normalized ion pressure and 
the adiabatic law P = N 3 , giving the initial condition 
T = N 2 . 

In Fig. 2, we present a comparison between the density 
profiles obtained from the fluid and Vlasov simulations, 
at different times. In the upper panel, the ion-electron 
temperature ratio 77 — 0.1, and the electron-ion density 
ratio a — 0.25. We see that both the fluid (left) and 
Vlasov (right) solutions exhibit shocks, where the shock 
front is distinct in the fluid solution and more diffuse in 
the Vlasov solution. The corresponding ion distribution 
function is displayed in Fig. 3. We observe that the for- 
mation of the shock at t = 3400 is located at z w 7000. 
It is associated with a "kink" in the distribution func- 
tion. A population of ions have also been accelerated by 
the shock. The middle panels of Fig 2 are for r\ = 0.1 
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FIG. 2: The profile of the ion density obtained from nu- 
merical solutions of the fluid equations (left panels) and the 
Vlasov equation (right panels). In the upper panels, the pro- 
files are shown at t = (dash-dotted lines), t = 1700 (dashed 
lines) and t = 3400 (solid lines) , while in the middle and lower 
panels the profiles are shown at t = 0, t = 500 and t = 1000. 
Parameters are r\ — 0.1 and a = 0.25 (upper panels), rj = 1 
and a = 0.25 (middle panels) and r\ = 1 and a = 0.05 (lower 
panels) . 
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FIG. 3: The ion distribution function at t = (upper panel), 
t = 1700 (middle panel) and t — 3400 (lower panel) as a 
function of x and v. Parameters are rj = 0.1 and a = 0.25. 



and a — 0.25. Here, the fluid solution exhibits clear 
shocks, while the Vlasov simulation shows only a phase 
of self-steepening at t = 500, followed by an expansion of 
the diffuse shockfont at t — 1000. The ion distribution 
function in Fig. 4 shows that the shockfront is strongly 
Landau damped for this case. Finally, the bottom panels 
of Fig. 2 show results for a = 0.05 and r] = 1. In this 
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FIG. 4: The ion distribution function at t — (upper panel), 
t = 500 (middle panel) and t = 1000 (lower panel) as a func- 
tion of x and v. Parameters are rj = 1 and a = 0.25. 
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FIG. 5: The ion distribution function at t — (upper panel), 
t = 500 (middle panel) and t = 1000 (lower panel) as a func- 
tion of x and t>. Parameters are rj — 1 and a = 0.05. 

case, the fluid solution again shows a shock in the front 
end of the pulse, but also the rear end of the shock steep- 
ens, which can be seen at z » —6000 for t = 1000. The 
steepening of the pulse for low-amplitude density per- 
turbations in the rear of the pules can be explained by 
that the wave speed decreases for small-amplitude den- 
sity perturbations (N < 1.1), as seen in Fig. 1, while 
it increases again for large-amplitude density perturba- 
tions. The Vlasov solution again shows a diffusive shock 
in the front, while it reproduces the steepening of the den- 
sity in the rear of the pulse. In the middle panel of Fig. 
5, the ion distribution function shows the self-steepening 
phase of the shockfront, and the lower panel shows the 
diffusion of the shock by shock- accelerated ions. In the 



rear end of the pulse, the distribution function forms a 
"kink," clearly seen at z » —6000 in the bottom panel. 

We have also performed simulations with smaller am- 
plitudes of the pulses (not shown here) and they exhibit 
essentially the same behavior as in the large-amplitude 
case. It is interesting to note that it is the strongly 
heated and shock-accelerated ions in the pulse that leads 
to Landau damping by overtaking the pulse. The heat- 
ing of the ions is due to adiabatic compression, lead- 
ing to a higher thermal speed of the ions inside the 
pulse than in the equilibrium plasma. Another effect 
is that the fluid (mean) velocity of the ions further ac- 
celerates the ions. For Landau damping to be unimpor- 
tant, we thus have the condition that the wave speed 
must be much larger than the sum of the ion thermal 
and fluid velocities. Inserting Eqs. (15) and (16) into 
the inequality A+ 3> Vt + u, where the scaled ion ther- 
mal speed Vt = (rjT) 1 ^ 2 ps r] 1 ^ 2 N, we obtain the con- 
dition N^ 2 [1/(N + a - 1) + 3777V] 1 / 2 > (rjN 2 ) 1 / 2 , or 
(l/[(iV+a-l)iV77]+3) 1 / 2 > 1. This condition is fulfilled 
if rj <C 1 (leading to the sharp shock seen in Fig. 3 and 
the upper right panel of Fig 2) or/ and if the electrons are 
evacuated due to the dust so that and at the same 

time N w 1. The latter corresponds to the case where 
the sign of the coefficient in front of the low-amplitude 
nonlinear term becomes negative, so that there will be 
a shock in the rear end of the pulse while the front of 
the shock expands, in agreement with the observations 
in Fig. 5 and lower right panel of Fig 2. The expansion 
of the shockfront at high dust densities has also been 
observed in the experiment 

To summarize, we have presented the dynamics of fully 
nonlinear, nondispersive dust ion acoustic waves in an 
unmagnetized dusty plasma. By using the Boltzmann 
electron distribution as well as the hydrodynamic equa- 
tions for the warm ion fluid and quasi-neutrality con- 
dition, we have represented the governing equations in 
the form of a master equation whose characteristics have 
been found analytically. The fluid equations has been 
solved to obtain the density and velocity profiles of the 
DIA shock waves, which exhibit the steepening of the 
waveforms both in the front and rear depending upon 
the values of a. We have also compared our theoretical 
results with those obtained from computer simulations of 
the time dependent Vlasov equation. The Vlasov solu- 
tion shows a diffuse shock in the front end of the pulse, 
due to strong Landau damping, while a sharp shock de- 
velops in the rear end of the pulse, similar to the results 
from the simulation of Eqs. (5) and (6). 
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